6  Linear Elasticity

Characteristic for continuous solids is that a solid can possess a non-trivial stress distribution at rest. The question often is, how such an arbitrary stress distribution translates into deformation and displacement of the material. Elasticity theory for continuum materials provides an explanation for a particular class of solids. It is a very old subject. The work of Hooke dates back to 1678. Elasticity theory has been developed in parallel to fluid mechanics, however, shares many analogies in the way the continuous material is described.

6.1 Hooke’s law

Assume, we have a rod of length \(l\), diameter \(s\), and cross-sectional area \(S\) that is subject to a load \(P\) in longitudinal direction. As a consequence the rod is strained, which implies a longitudinal stretching of \(\Delta l\) and a cross-sectional thinning of \(\Delta S\). Classical Hooke’s law states that the stretching can be determined from

\[ \underbrace{\frac{\Delta l}{l}}_{\text{strain } \epsilon } = \underbrace{\frac{1}{E}}_{\substack{\text{proportionality} \\ \text{constant}}} \; \underbrace{\frac{P}{S}}_{\text{stress } \sigma }. \]

Hence, internal stress is related to longitudinal (axial) strain within the rod:

\[ \sigma = E \epsilon. \]

In particular, Hooke’s law states that stress and strain are proportional with a proportionality constant given by the Young modulus \(E\) that characterizes a mechanical stiffness of a material. Accordingly, we have

\[ \frac{\Delta S}{S} = - \nu \frac{\Delta l}{l} = - \nu \frac{1}{E}\,\frac{P}{S}, \]

meaning that the transverse strain \(\tfrac{\Delta S}{S}\) in direction perpendicular to the loading differs from the longitudinal strain \(\tfrac{\Delta l}{l}\) by a factor denoted as the Poisson ratio \(\nu\).

So far, we looked at one particular rod and one specific loading. Hooke’s law now states that the two material parameters (Young modulus \(E\) and Poisson’s ratio \(\nu\)) are - within a clearly defined scale and load regime referred to as the linear elastic regime - is independent of the rod’s dimension \(l\) and \(S\) as well as of the value of the load \(P\).

In this lecture we will focus on situations, in which \(E\) is very large, e.g. \(\mathcal O ( \text{GPa})\)), such that longitudinal stretching is very small in comparison to the original length of the rod (\(\Delta l \ll l\)). This means that the strain is very small: \(\tfrac{\Delta l }{l} \ll 1\). For instance, if we load a steel rod of 1m length and 1cm in diameter with a mass of 1t in longitudinal direction, we will observe a longitudinal stretching in the order of \(\mathcal O(10^{-3})\) and a cross-sectional thinning in the order of \(\mathcal O(10^{-6})\). The steel rod retains its original configuration, once the load is released.

Such a material is referred to as a linear, elastic solid. The question now is: How do we carry over Hooke’s law to a continuous material body? Essentially, we are interested in the displacement of material particles in response to arbitrary, potentially time dependent loading. Can we re-write formerly derived balance laws, e.g. momentum balance \(\rho \tfrac{d}{dt} \mathbf v = \nabla \cdot \sigma\) to allow for the modeling of displacements within a solid material? How can we formulate a generalized version of Hooke’s law as a closure for the momentum balance’s right hand side \(\nabla \cdot \sigma\)? Answering these questions constitutes the roadmap for this section on linear elastic solids.

6.2 Deformation kinematics

Let’s assume, we have a material that is subject to a loading and hence deformed. Deformation is given by a mapping \[\phi : B \rightarrow B',\] also referred to as the deformation function which maps a materials original (reference) configuration \(\mathbf \zeta \in B\) is mapped onto a current (deformed) material configuration \(\mathbf x = \phi(\mathbf \zeta,t) \in B'\). Points in \(B\) are denoted by \(\mathbf \zeta = (\zeta_1,\zeta_2,\zeta_3)^T\) and refer to material coordinates, whereas points in \(B'\) are denoted by \(\mathbf x=(x_1,x_2,x_3)^T = (\phi_1(\mathbf \zeta,t),\phi_2(\mathbf \zeta,t),\phi_3(\mathbf \zeta,t))^T\) refer to deformed or spatial coordinates.

Deformation mapping \(\phi\) codes local changes in shape, for instance local stretching, which is referred to as strain. Consider an open ball \(\Omega \subset B\) of radius \(\epsilon\) centered around \(\mathbf \zeta_0\). \(\phi(\Omega)\) then refers to the ball after deformation

\[\phi(\Omega) = \{ \mathbf x \in B' | \mathbf x = \phi(\mathbf \zeta,t), \mathbf \zeta \in B \}.\]

Any relative difference of \(\Omega\) versus \(\phi(\Omega)\) as \(\epsilon \rightarrow 0\) is called strain at \(\mathbf \zeta_0\). Note, however, that \(B'\) not necessarily indicates a changes in shape. It can likewise code a rotation or a translation of reference configuration \(B\).

The type of deformation can be further characterized by the deformation gradient \(F\) defined as

\[ F := \nabla \phi(\mathbf \zeta) = \frac{\partial \phi}{\partial \mathbf \zeta} \qquad \Leftrightarrow \qquad F_{ij}:=\frac{\partial \phi_i}{\partial \zeta_j}. \]

It is hence given as the derivative of each component of the deformed configuration \(\mathbf x\) with respect to each component of the reference position vector \(\mathbf \zeta\). \(F\) provides information on the local behaviour of the deformation mapping. By means of a Taylor expansion, we get

\[\triangle \mathbf x = \mathbf x - \mathbf x_0 = F(\mathbf \zeta_0) \triangle \mathbf \zeta + \mathcal O(|\triangle \mathbf \zeta|^2).\]

We note the following classics:

Pure rigid body displacement yields \(F=I\).

Example: \(x_1 = \zeta_1 + 5, x_2 = \zeta_2, x_3 = \zeta_3.\)

Pure rigid body rotation yields a deformation gradient corresponding to a rotation matrix.

Example: \(x_1 = \cos \theta \zeta_1 - \sin \theta \zeta_2, x_2 = \sin \theta \zeta_1 + \cos \theta \zeta_2, x_3 = \zeta_3\).

Pure stretching in axial direction yields a diagonal deformation gradient \(F = diag(\alpha_1,\alpha_2,\alpha_3)\).

Example: \(x_1 = \alpha_1 \zeta_1, x_2 = \alpha_2 \zeta_2, x_3 = \alpha_3 \zeta_3\)

Pure shear yields a symmetric deformation gradient \(F^T = F\).

6.3 Polar decomposition of \(F\)

In order to decompose \(F\) systematically, we introduce the polar decomposition. Our aim is to systematically characterize the deformation mapping \(\phi\) based on our knowledge of the deformation gradient \(F\).

\(F\) can be uniquely decomposed into

\[ F = Q \cdot U \qquad \text{or alternatively} \qquad F = V \cdot Q. \]

Here, \(Q\) is an orthogonal matrix and \(U\) and \(V\) are symmetric, positive definite matrices respectively. If we interpret \(F\) as a linear operator this means that \(F\) is decomposed into subsequently applied operators \(U\) and \(Q\), or \(Q\) and \(V\).

Proof: We define

\[ U := (F^T F)^{\tfrac{1}{2}} \qquad \text{and} \qquad Q:= F U^{-1}, \]

and deduce the following:

  1. \(F^T F\) is symmetric
  2. \(\mathbf x^T F^T F \mathbf x = (F \mathbf x)^T F \mathbf x > 0\) for all \(\mathbf x\).
  3. Hence, \(F^T F\) is spd and has a full set of eigenvalues. It can be diagonalized and square rooted.
  4. As a consequence, \(U\) exists and is invertible, such that also \(Q\) exists.
  5. Finally, we can show that \(Q\) is orthogonal (\(Q^T Q = I\)), which completes the proof.

The same can be shown for

\[ V := (F F^T)^{\tfrac{1}{2}} \qquad \text{and} \qquad Q:= V^{-1} F, \]

\(U\) and \(V\) appearing in the right and left polar decomposition of \(F\) code information on the strain and are referred to as the right and left stretch tensors. They are not identical in general, but share the same eigenvalues, referred to as the principle stretches. The sets of eigenvectors \(\mathbf u_i\) and \(\mathbf v_i\) are called the right and left principle stretch directions.

The deformation gradient’s polar decomposition can be physically interpreted as subsequent rotation \(Q\) and stretch \(U\) and \(V\) or subsequent stretch and rotation respectively. We get the following transformation rules for a material particle \(\Delta \mathbf \zeta\):

\[ \begin{aligned} \Delta \mathbf x &= F \Delta \mathbf \zeta = Q ( U \mathbf \Delta \zeta)\\ \Delta \mathbf x &= F \Delta \mathbf \zeta = V ( Q \mathbf \Delta \zeta). \end{aligned} \]

Remark 2

While the stretch tensors differ in general (\(U\neq V\) in general), we have the same rotation \(Q\) in both cases. Strain before and after rotation are related:

\[V=Q U Q^T.\]

There are other relevant strain tensors. For example:

  • \(B : = V^2 = F F^T\): left Cauchy-Green tensor

  • \(C : = U^2 = F^T F\): right Cauchy-Green tensor

Note, that these have been used for the construction of \(U\) and \(V\).

Rigid body deformation means vanishing strain: \[U=V=B=C=I.\]

Volume preserving deformation is characterized by \(\text{det}(F) = 1\), because

\[ \int_{\Omega_{\mathbf x}} d \mathbf x = \int_{\Omega_{\mathbf \zeta}} \underbrace{\left| \tfrac{\partial \mathbf x}{\partial \mathbf \zeta}\right|}_{\mathclap{\substack{{\text{det}(F) \,:\, \text{Jacobi}} \\ \text{determinant of the} \\ \text{deformation gradient}}}} d \mathbf \zeta. \]

With that, we can single out isotropic dilation and volume preserving strain and yield:

\[ F= \underbrace{Q}_{\mathclap{\text{rotation}}} \cdot \underbrace{\text{det}(F) I }_{\substack{\mathclap{\text{isotropic}} \\ \mathclap{\text{dilation}}}} \cdot \underbrace{\dfrac{1}{\text{det}(F)} U}_{\substack{\mathclap{\text{volume}} \\ \mathclap{\text{preserving strain}}}}. \]

This means that locally any deformation can be understood as subsequent volume preserving stretching, isotropic dilation and rotation.

6.4 General material theory in a nutshell

In elasticity theory we are interested in materials, whose stress state depends on (local) strain alone. This can be further formalized as follows:

6.4.1 Elasticity

We call a solid material elastic, if its constitutive law \(\sigma = \hat \sigma (F)\) depends only on the deformation gradient \(F\).

6.4.2 Symmetry

If we transform the material’s reference configuration according to

\[\mathbf \zeta \to \mathbf \eta ( \mathbf \zeta ),\]

we get

\[ F_{ij}^{(\eta)} = \tfrac{\partial \phi_i}{\partial \eta_j} = \tfrac{\partial \phi_i}{\partial \zeta_k} \tfrac{\partial \zeta_k}{\partial \eta_j} = F_{ik}^{(\zeta)} \underbrace{P_{kj}}_{\mathclap{\text{transformation matrix}}} \qquad \Leftrightarrow \qquad F{(\eta)} = F{(\zeta)} \cdot P \]

The material is referred to as symmetric with respect to \(P\), if the result of the constitutive relation \(\hat \sigma (*)\) in unaffected by the change of the reference configuration, hence

\[ \hat \sigma^{(\zeta)} (*) = \hat \sigma^{(\eta)} (* P). \]

6.4.3 Isotropy

A material that is symmetric in the above sense for all \(P=Q\) with \(Q\) orthogonal is called isotropic. Note, that if we combine this with the former definition of elasticity, isotropy implies

\[ \hat \sigma (F) \underbrace{=}_{\text{polar decomposition}} \hat \sigma ( V \cdot Q) \underbrace{=}_{\text{symmetry}} \hat \sigma ( V \cdot Q \cdot P) \underbrace{=}_{\text{isotropy: pick } P = Q^T} \hat \sigma ( V). \]

This means, that in an isotropic, elastic material, the constitutive relation \(\hat \sigma (*)\) is a function of rotation-less strain alone!

6.4.4 Objectivity

Objectivity, finally, states that the constitutive relation is invariant under Euclidian transformations!

In order to understand this, let’s consider a Euclidian transformation of the deformed configuration

\[ \mathbf y = Q(t) \mathbf x + \mathbf b(t), \]

in which \(Q(t)\) is orthogonal and denotes a (potentially time-dependent) rotation, and \(\mathbf b(t)\) denotes a (potentially time dependent) shift. Similarly to before, we find that

\[ \begin{aligned} F^{(\mathbf y)} = Q \cdot F^{(\mathbf x)} \qquad \text{and hence } \qquad \sigma^{(\mathbf y)} = Q \cdot \sigma^{(\mathbf x)} \cdot Q^T. \end{aligned} \]

In general, this would imply that also the constitutive relation has to be transformed, namely $\(\hat \sigma (F^{(\mathbf x)}) \neq \hat \sigma^{(\mathbf y)} (F^{(\mathbf y)}) = \hat \sigma (Q \cdot F^{(\mathbf x)})\).

Objectivity, however, yields

\[ \hat \sigma^{(\mathbf y)} = \hat \sigma^{(\mathbf x)}. \]

6.5 Displacement kinematics

After a load has been applied, a material particle at position \(\mathbf \zeta\) in reference configuration \(B\) is displaced according to

\[ \mathbf w ( \mathbf \zeta) := \mathbf x - \mathbf \zeta = \phi ( \mathbf \zeta,t) - \mathbf \zeta,\]

which is referred to as the displacement field. Partial differentiation of the displacement field with respect to the material coordinate \(\mathbf \zeta\) yields the (material) displacement gradient

\[ \underbrace{H}_{\substack{\mathclap{\text{displacement}} \\ \mathclap{\text{gradient}}}} : = \nabla_{\mathbf \zeta} \phi(\mathbf x) - I = \tfrac{\partial \phi_i}{\partial \zeta_j} - \delta_{ij} = \underbrace{F}_{\substack{\mathclap{\text{deformation}} \\ \mathclap{\text{gradient}}}} - I \]

We can also solve this for the deformation gradient, hence \(F = H + I\).

Now recall our previous result regarding the deformation gradient’s polar decomposition, when we singled out strain from rotations! In order to construct the strains, e.g. \(V\), we first had to compute Cauchy-Green tensors. Let’s look at the left Cauchy-Green tensor:

\[ B_{ij} = F_{ik} F_{kj} = \left( \delta_{ik} + \frac{\partial}{\partial \zeta_k} w_i \right) \left( \delta_{kj} + \frac{\partial }{\partial \zeta_j} w_k \right) = \underbrace{\delta_{ij} + \frac{\partial }{\partial \zeta_j} w_i + \frac{\partial}{\partial \zeta_i}w_j}_{\text{linear contribution: } I + H + H^T} + \underbrace{\frac{\partial }{\zeta_k \zeta_k} w_i w_j}_{\substack{\mathclap{\text{non-linear}} \\ \mathclap{\text{contribution}}}}. \]

This means that the Cauchy-Green tensor consists of a constant part, a part that is linear in the displacement gradient and further non-linear contributions:

\[ B = F \cdot F^T = I + H + H^T + \text{non-linear contributions}. \]

This fact explains the complicated (because typically non-linear) nature of general solid elasticity models. Here, however, we want to follow the idea of Hooke’s law focus on material behaviour that is characterized by very small displacements. This simplifies the situation significantly. as discussed in the next section.

6.6 Small displacement fields

We assume that the displacement variation within the material is very small and formalize this via \(\tfrac{\partial w_i}{\partial \zeta_j} = \mathcal O(\epsilon)\) with \(\epsilon \ll 1\). As a consequence to displacements varying only slightly, we find that the deformation gradient \(F\) varies only slightly from identity
\[ F = I + \mathcal O(\epsilon). \]

Let’s furthermore assume that we want to keep all terms that scale with \(\mathcal O(\epsilon)\), yet want to neglect terms that scale with \(\mathcal O(\epsilon^2)\). The latter assumption is referred to as the small displacement assumption. The small displacement assumption leads to a significant simplifications:

First, we find that up to order \(\mathcal O(\epsilon)\) there is no difference between a differentiation of the displacement field with respect to coordinates of the reference frame \(x_i\) or material coordinates (\(\zeta_i\)):

\[ \begin{aligned} \frac{\partial w_i}{\partial \zeta_j} = \frac{\partial w_i}{\partial \phi_k} \frac{\partial \phi_k}{\partial \zeta_j} = \frac{\partial w_i}{\partial x_k} F_{kj} = \frac{\partial w_i}{\partial x_k} \left( \delta_{kj} + \frac{\partial w_k}{\partial \zeta_j}\right) = \frac{\partial w_i}{\partial x_j} + \mathcal O(\epsilon^2). \end{aligned} \]

This means basically, that the difference between a Eulerian and Lagrangian description of displacement field is of order \(\mathcal O(\epsilon^2)\) and can be neglected.

We furthermore find that left and right Cauchy-Green tensor are identical with vanishing non-linear contributions \[ \begin{aligned} B = C = I + H + H^T + \mathcal O(\epsilon^2), \end{aligned} \]

which implies that \(U\) and \(V\) are also identical given by \[ \begin{aligned} U = V = I + \frac{1}{2} (H + H^T) + \mathcal O(\epsilon^2), \end{aligned} \]

as we get \(U^T U = B_{ij}\).

The inverse stretch tensors accordingly simplify into \[ \begin{aligned} U^{-1} = V^{-1} = I - \frac{1}{2} \left( H + H^T \right) + \mathcal O(\epsilon^2). \end{aligned} \]

Finally, the orthogonal rotation tensor can be determined to be \[ \begin{aligned} Q = I + \frac{1}{2} \left( H - H^T \right) + \mathcal O(\epsilon^2). \end{aligned} \]

We conclude that for small displacements, the polar decomposition is given as \[ \begin{aligned} F = \underbrace{I + \frac{1}{2} \left( H + H^T \right)}_{\substack{\mathclap{\text{corresponds to}} \\ \mathclap{\text{strain } V}}} + \underbrace{\frac{1}{2} \left(H - H^T \right)}_{\substack{\mathclap{\text{corresponds to}} \\ \mathclap{\text{rotation } Q}}} + \mathcal O(\epsilon^2). \end{aligned} \]

The symmetric tensor \(\frac{1}{2} \left( H + H^T \right)\) is also referred to as the small strain tensor. We find that

  1. Rigid body deformation yields \(B=C=I\), hence \(\tfrac{1}{2} \left( H + H^T \right)=0\)

  2. Rotation is controlled by \(\tfrac{1}{2} \left( H - H^T \right)=0\)

  3. Volume preservation implies \(\text{det}(F) = 1\), hence \(\tfrac{\partial w_k}{\partial x_k}=0\)

which reminds us of findings, when we decomposed the velocity vector (see Section 4.2).

6.7 Continuous formulation of Hooke’s law

In the small displacement limit, the constitutive relation for a linearly elastic, isotropic, objective solid is given by Hooke’s law:

\[ \sigma = \lambda tr (D) \mathbf I + 2 \mu D, \]

with small-strain tensor \(D\) being given the symmetric part of the displacement gradient \(D:=\tfrac{1}{2}(H + H^T)\). Parameters \(\lambda\) and \(\mu\) denote the so-called Lamé constants and can be measured for each material. Note, that \(tr(D) = \partial_k w_k\) codes the material’s change in volume.

Remark

The sign convention differs between classical fluid dynamics and classical continuum mechanics for solids. Here, we chose the solid convention, assuming compressive stresses to be negative.

Hooke’s law can alternatively be stated in terms of the small-strain deviator

\[D^{(d)} := D - \frac{1}{3} tr (D)I,\]

for which \(tr (D^{(d)}) = 0\). This yields

\[ \begin{aligned} \sigma &= 2 \mu D^{(d)} + K tr (D) I, \end{aligned} \]

with \(K= \lambda + \tfrac{2}{3} \mu\). In index notation, this reads

\[ \begin{aligned} \sigma_{ij} &= 2 \mu \left( D_{ij} - \frac{1}{3} D_{kk} \delta_{ij} \right) + K D_{kk} \delta_{ij}. \end{aligned} \tag{6.1}\]

The trace can be computed via convolution (setting \(i=j\) and summation), which results in

\[ \sigma_{kk} = 3 K D_{kk} \qquad \Leftrightarrow \qquad tr(\sigma) = 3 K tr(D). \tag{6.2}\]

Hence, the trace of the stress tensor and the trace of the small-strain tensor are proportional to each other.

Let’s for instance consider a hydrostatic situation (lithostatic/ glaciostatic), in which \(\sigma = - p I\), hence \(tr(\sigma) = - 3 p\) and consequently

\[ tr(D) = - \frac{p}{K}. \]

That means that the volume strain \(tr(D)\) is proportional to pressure p with proportionality constant \(-K^{-1}\). \(K\) is therefore called the volume compression modulus.

For shear stresses, we get \(\delta_{ij} = 0\), hence

\[ \sigma_{ij} = 2 \mu D_{ij} = \mu ( \partial_j w_i + \partial_i w_j). \]

Parameter \(\mu\) is hence referred to as the shear modulus.

At the beginning of this chapter, we started from classical Hooke’s law written in terms of Youngs modulus and Poisson’s ratio. The question now is: How are both versions of Hooke’s law connected?

In order to analyze this further, we solve Equation 6.1 for the small-strain \(D\) by first substituting in Equation 6.2

\[ \sigma_{ij} = 2 \mu \left( D_{ij} - \frac{1}{9K} \sigma_{kk} \delta_{ij} \right) + \frac{1}{3} \sigma_{kk} \delta_{ij}. \]

Now, it is possible to solve for \(D_{ij}\), which yields

\[ D_{ij} = \frac{1}{9K} \sigma_{kk} \delta_{ij} + \frac{1}{2 \mu} ( \sigma_{ij} - \frac{1}{3} \sigma_{kk} \delta_{ij} ). \]

In our initial strained rod example, we had \(D_{11} = \tfrac{\Delta l}{l}\) and \(D_{22} = \tfrac{\Delta S}{S}\). Axial load was exerted in \(x_1\) direction, hence \(\sigma_{11} = \tfrac{P}{S}\). It was zero elsewhere. This resolves in

\[ \begin{aligned} D_{11} = \frac{\Delta l}{l} &= \frac{1}{9K} \frac{P}{S} + \frac{1}{2 \mu} \left( \frac{P}{S} -\frac{1}{3}\frac{P}{S} \right) \\ &= \left( \frac{1}{9K} + \frac{1}{3 \mu} \right) \frac{P}{S} \overset{!}{=} \frac{1}{E}\frac{P}{S} \end{aligned} \] \[ \begin{aligned} D_{22} = \frac{\Delta S}{S} &= \frac{1}{9K} \frac{P}{S} - \frac{1}{6 \mu} \frac{P}{S} \overset{!}{=} - \frac{\nu}{E}\frac{P}{S} \end{aligned} \]

We can solve these relations for \(\mu, K\), and \(\lambda\) and find

\[ \mu=\frac{E}{2(1+\nu)} \qquad K=\frac{E}{3(1-2\nu)} \qquad \lambda=\frac{\nu E}{(1+\nu)(1-2\nu)}. \]

It provides us with a relation between the Lamé constants and Young’s modulus, and Poisson’s ratio respectively.

6.8 Theory of linear elasticity

The major goal of this chapter is to derive a mathematical process model, hence a PDE system, that allows us to compute the spatio-temporal evolution of the displacement field \(\mathbf w\).

We start from the momentum balance (body forces are neglecting for time being):

\[ \rho \frac{D}{Dt} \mathbf v = \nabla \cdot \sigma \]

First, we rewrite total derivative in Lagrangian coordinates

\[\frac{D}{Dt} \mathbf v = \frac{\partial}{\partial t} \mathbf{\tilde v},\]

which allows us to substitute in the deformation mapping and displacement field:

\[ \frac{\partial}{\partial t} \underbrace{\mathbf{\tilde v}(\zeta,t)}_{\substack{\mathclap{\text{velocity in}} \\ \mathclap{\text{Lagrangian reference}}}} = \frac{\partial}{\partial t} \left(\frac{\partial}{\partial t} \phi(\zeta,t)\right) = \frac{\partial^2}{\partial t^2} \left( \mathbf{w}(\zeta,t) + \zeta \right) = \frac{\partial^2}{\partial t^2} \mathbf{w}(\zeta,t). \]

This finally yields

\[ \rho \frac{\partial^2}{\partial t^2} \mathbf w = \nabla \cdot \sigma. \]

Remark

Recall, that we found earlier that in the small displacement limit, Eulerian and Lagrangian spatial derivatives of the displacement tensor are identical up to order \(\epsilon^2\). This means that we do not have to differentiate between both reference frames when evaluating the right hand side \(\nabla \cdot \sigma\).

We can now substitute in Hooke’s law in continuous materials: \[ \rho \frac{\partial^2}{\partial t^2} \mathbf w = 2 \mu \nabla \cdot D + \lambda \nabla \cdot (tr(D) I). \]

in terms of small strain tensor \(D = \tfrac{1}{2} \left( H + H^T\right)\). It’s divergence is given by

\[ \begin{aligned} \nabla \cdot D &= \frac{1}{2} \nabla \cdot \left( H + H^T \right) \\ &= \frac{1}{2} \nabla \cdot \left( \nabla \mathbf w + (\nabla \mathbf w)^T \right) = \frac{1}{2} \left( \nabla (\nabla \cdot \mathbf w) + \Delta \mathbf w \right). \end{aligned} \]

Exploiting the operator identity

\[ \nabla \cdot ( tr(D) I) = \nabla (\nabla \cdot \mathbf w) \]

finally gets us the dynamic theory of linear elasticity

\[ \rho \frac{\partial^2}{\partial t^2} \mathbf w = \mu \Delta \mathbf w + \frac{\mu}{1-2 \nu} \nabla (\nabla \cdot \mathbf w). \]

Remark

Note that factor \(\frac{\mu}{1-2 \nu}\) is found by substituting in the previously derived relations for \(\lambda\)!

With the vector identity \(\nabla (\nabla \cdot \mathbf w) = \Delta \mathbf w + \nabla \times \nabla \times \mathbf w\), this can alternatively be written as

\[ \rho \frac{\partial^2}{\partial t^2} \mathbf w = \frac{\mu}{1-2\nu} \left( 2 ( 1-\nu) \Delta \mathbf w + \nabla \times \nabla \times \mathbf w \right). \]

6.9 Static theory of elasticity

The static theory of elasticity goes back to Navier (1827). It results from the dynamic theory by neglecting transient effects, hence

\[\frac{\partial^2}{\partial t^2} \mathbf w = 0,\]

which yields

\[ 0 = ( 1- 2 \nu) \Delta \mathbf w + \nabla (\nabla \cdot \mathbf w). \]

The static theory has important properties:

  1. The divergence of the displacement field \(\nabla \cdot \mathbf w\), namely \(tr(D)\), is a harmonic function, hence a solution to the homogeneous Laplace equation \(\Delta ( \cdot ) = 0\). We see this from \[ \begin{aligned} 0 &= \nabla \cdot \left( ( 1- 2 \nu) \Delta \mathbf w + \nabla (\nabla \cdot \mathbf w) \right) \\[1em] &= \underbrace{( 1-2\nu) \Delta \nabla \cdot \mathbf w}_{\nabla \cdot \text{ and } \Delta \text{ commute}} + \underbrace{\nabla \cdot \nabla (\nabla \cdot \mathbf w) }_{= \Delta \nabla \cdot \mathbf w} \\[1em] &= 2( 1-\nu) \Delta (\nabla \cdot \mathbf w) \end{aligned} \]

  2. The displacement field \(\mathbf w\) is biharmonic, which can be seen by applying the Laplacian operator to static elasticity theory: \[ \begin{aligned} 0 &= \Delta \left( ( 1- 2 \nu) \Delta \mathbf w + \nabla \nabla \cdot \mathbf w \right) \\[1em] &= (1-2\nu) \Delta \Delta \mathbf w + \underbrace{\Delta \nabla \nabla \cdot \mathbf w }_{=\nabla \underbrace{\Delta \nabla \cdot \mathbf w}_{=0} =0} \\ &= ( 1-2\nu) \Delta \Delta \mathbf w, \end{aligned} \] hence \(\Delta \Delta \mathbf w = 0\). The latter defines a biharmonic function.

This offers a pathway towards a numeric solution: In order to solve for the displacement field, we finally need boundary conditions, e.g.

  • a prescribed displacement given as \(\mathbf w = f(\mathbf x)\) for \(\mathbf x \in \partial \Omega\), or

  • a prescribed traction vector \(\sigma \cdot \mathbf n = f(\mathbf x)\) for \(\mathbf x \in \partial \Omega\).

6.10 Wave propagation

Wave propagation is a special case of dynamic elasticity. Let’s decompose the displacement field into a divergence free part referred to as solenoidal part, and a curl free part referred to as the irrotational part. This type of decomposition is very famous and known as the Helmholtz decomposition: \[ \mathbf w ( \mathbf x,t) = \underbrace{\mathbf w_1(\mathbf x,t)}_{\nabla \cdot \mathbf w_1=0} + \underbrace{\mathbf w_2(\mathbf x,t)}_{\nabla \times \mathbf w_2=0} \tag{6.3}\]

Substitution into the dynamic wave equation yields

\[ \rho \frac{\partial^2}{\partial t^2} (\mathbf w_1 + \mathbf w_2) = \mu \Delta (\mathbf w_1 + \mathbf w_2) + \frac{\mu}{1 - 2 \nu} \nabla (\nabla \cdot (\mathbf w_1 + \mathbf w_2)). \]

Applying the curl operator to the latter and exploiting vector identities yields

\[ \nabla \times \underbrace{\left( \rho \frac{\partial^2}{\partial t^2} \mathbf w_1 - \mu \Delta \mathbf w_1 \right)}_{=(I)} = 0. \]

Note, that per construction of the Helmholtz decomposition also the divergence of the resulting expression vanishes: \(\nabla \cdot (I) = 0\).

We now apply the divergence operator to the decomposed wave equation and find

\[ \nabla \cdot \underbrace{\left( \rho \frac{\partial^2}{\partial t^2} \mathbf w_2 - \frac{2(1-\nu)}{1-2\nu}\mu \Delta \mathbf w_2 \right)}_{=(II)} = 0, \]

for which we similarly get \(\nabla \times (II) = 0\).

Both \((I)\) and \((II)\) have a vanishing divergence and a vanishing curl. As a consequence, both have to posses a scalar potential, which is a harmonic function \[ \begin{aligned} (I) = \nabla \phi_I \quad &\text{with}& \Delta \phi_I = 0 \\ (II) = \nabla \phi_{II} \quad &\text{with}& \Delta \phi_{II} = 0. \end{aligned} \]

For a body in which displacements vanish in the far-field, we hence conclude \(\phi_{I} = \phi_{II}=0\) to be the only solution. This in return means that both \(\mathbf w_1\), and \(\mathbf w_2\) fulfill the wave equations:

\[ \begin{aligned} \frac{\partial^2}{\partial t^2} \mathbf w_1 - a_1 \Delta \mathbf w_1 = 0 \quad \text{where} \quad &a_1 = \frac{\mu}{\rho} \\ \frac{\partial^2}{\partial t^2} \mathbf w_2 - a_2 \Delta \mathbf w_2 = 0 \quad \text{where} \quad &a_2 = \frac{2(1-\nu)}{1-2\nu} \frac{\mu}{\rho} > a_1. \end{aligned} \]

Waves of the first type have no change in volume (\(\nabla \cdot \mathbf w_1=0\)). They are characterized by distorsion of particles. Alternative names are shear waves, transverse waves, solenoidal waves or S-waves.

Waves of the second type produce volume changes and are referred to as compression expansion waves, pressure waves, longitudinal, dilational, irrotational waves or P-waves.